function [output] = CRS_main(input_array)
% This function uses a sequence of clustering structure to perform 
% Canay, Romano, and Shaikh (2017). It produces a vector of test restuls of 
% 0 or 1. "testCRS_star" is the test result using the clustering selected 
% according to a size-power tradeoff criterion. 


%% data input
method = input_array.method;
data = input_array.data;
Sigma_hat = input_array.Sigma_hat;
group_matrix = input_array.group_matrix_km;
alpha_sig = input_array.alpha_sig;
alt_vec = input_array.alt_vec;
alt_vec_tradeoff = input_array.alt_vec_tradeoff;


%% adjust p-value thresholds to control sizes for each clustering
G_vec = max(group_matrix); % vector of group numbers
l_g = length(G_vec);
% t1_vec = zeros(1,length(G_vec));
ave_power_vec = zeros(1,length(G_vec));
p_adj_vec = zeros(1,length(G_vec));
rej_rate_vec = zeros(1,length(G_vec));
for jj = 1:numel(G_vec)
    
    membership = group_matrix(:,jj);
    
    % find p-value threshold such that rejection rate is nominal
    [rej_rate,p_adj] = rej_rate_CRS_mean(Sigma_hat,membership,...
        alpha_sig,0);
    rej_rate_vec(jj) = min(rej_rate,alpha_sig);
    p_adj = min(p_adj,alpha_sig); % cannot be more aggressive
   
    power_vec = zeros(length(alt_vec_tradeoff),1);
    for i_alt = 1 : length(alt_vec_tradeoff)
        power_vec(i_alt) =  rej_rate_CRS_mean(Sigma_hat,...
            membership,p_adj,alt_vec_tradeoff(i_alt));
    end
    ave_power_vec(jj) = mean(power_vec);
    
    p_adj_vec(jj) = p_adj;
    
end
[ave_power,ind_G] = max(ave_power_vec);
group_star = group_matrix(:,ind_G);
p_adj = p_adj_vec(ind_G);
G_star = G_vec(ind_G);
estimated_size = rej_rate_vec(ind_G);


%% ESTIMATION
D = data.D;
X = data.X;
Y = data.Y;
        
switch lower(method)
    case 'ols'
        Z = D;
        
        alt_power_curve = linspace(alt_vec(1),alt_vec(end),50);
        
    case 'iv'
        Z = data.Z;
        
        alt_power_curve = linspace(alt_vec(2),alt_vec(end-1),50);
end

[bfmtmp,~,beta_g] = FamaMacbeth(D,X,Y,Z,group_star); % FM results
thetaHat_fm = bfmtmp(1);


%% P-VALUE
beta0_g = beta_g(:,1);

% generate random transformations for re-use
m_dBin = 1000; % number of simulated transformations
n_dBin = max(G_vec); % number of groups
d = ((rand(m_dBin,n_dBin)>.5)-.5)*2; % simulated transformations

[~,p_null] = CRS(beta0_g,d(:,1:max(group_star)),p_adj);

p_alt = 0*alt_vec;
for i_alt = 1 : length(alt_vec)
    beta0_g_alt = beta0_g+alt_vec(i_alt);
    [~,p_alt_i] = CRS(beta0_g_alt,d(:,1:max(group_star)),p_adj);
    p_alt(i_alt) = p_alt_i;
end


%% OUTPUT for estimation and inference table
output.estimate = thetaHat_fm;
output.p_adj = p_adj;
output.p_null = p_null;
output.p_alt = p_alt;


%% OUTPUT for power curve
% alt_power_curve = linspace(alt_vec(1),alt_vec(end),50);
p_power_curve = 0*alt_power_curve;
for i_alt = 1 : length(alt_power_curve)
    beta0_g_alt = beta0_g+alt_power_curve(i_alt);
    [~,p_alt_i] = CRS(beta0_g_alt,d(:,1:max(group_star)),p_adj);
    p_power_curve(i_alt) = p_alt_i;
end

output.alt_power_curve = alt_power_curve;
output.p_power_curve = p_power_curve;


%% p-value with fixed G

p_fixed_G = zeros(1,length(G_vec));
for i_g = 1 : l_g

    group = group_matrix(:,i_g);

    [~,~,beta_g] = FamaMacbeth(D,X,Y,Z,group); % FM results
    beta0_g_temp = beta_g(:,1); % group-level estimates
    [~,p_value] = CRS(beta0_g_temp,d(:,1:max(group)),...
        p_adj_vec(i_g)); % p-value
    p_fixed_G(i_g) = p_value;

end


%% OUTPUT for clustering table
output.G_star = G_star;
output.ave_power = ave_power;
output.p_fixed_G = p_fixed_G;
output.p_threshold_fixed_G = p_adj_vec;
output.estimated_size = estimated_size;

p_tradeoff = 0*alt_vec_tradeoff;
for i_alt = 1 : length(alt_vec_tradeoff)
    beta0_g_alt = beta0_g+alt_vec_tradeoff(i_alt);
    [~,p_alt_i] = CRS(beta0_g_alt,d(:,1:max(group_star)),p_adj);
    p_tradeoff(i_alt) = p_alt_i;
end
output.p_tradeoff = p_tradeoff;




